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Abstract 

A method of solving the time-dependent Schrodinger equation is pre- 
sented, in which a finite region of space is treated explicitly, with the bound- 
ary conditions for matching the wave-functions on to the rest of space re- 
placed by an embedding term added on to the Hamiltonian. This time- 
dependent embedding term is derived from the Fourier transform of the 
energy-dependent embedding potential, which embeds the time-independent 
Schrodinger equation. Results are presented for a one-dimensional model of 
an atom in a time-varying electric field, the surface excitation of this model 
atom at a jellium surface in an external electric field, and the surface excita- 
tion of a bulk state. 



1 Introduction 

The availability of ultra-short laser pulses [1] opens up new ways of studying 
time-dependent electronic processes in atoms fl2l|3l, molecules [4J and at solid 
surfaces 0. In surface physics the time delay of a core electron emitted through 
the surface potential barrier can now be compared on an attosecond timescale 
with valence electron photoemission in time-resolved photoemission experiments 
0. This makes it important to develop appropriate theoretical tools for studying 
such processes, in particular for solving the time-dependent Schrodinger equation 
accurately. A problem arises with the boundary in solving this equation - if an 
electron is ejected from an atom, for example, how is its wave-function treated as it 
propagates outwards towards the edge of the computational region? One approach 
is to apply absorbing boundary conditions at the edge of the region [7], but this 
is an approximation [8]; complex coordinate methods can also be used to remove 
the outgoing wave in solutions of the time-dependent Schrodinger equation [0. 
Other methods are used in atomic physics to study photoionization in particular - 
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the wave-function can be expanded in complex basis functions iTTOl . and Floquet 
methods use a Floquet-Fourier series expansion of the wave-function [fTTTl . 

In recent years there has been work on transparent boundary conditions, so 
that the solution of the time-dependent Schrodinger equation in some restricted 
region of interest, which we call region I, propagates out through the boundary 
into the rest of the system, region II, without reflection. It has been shown by 
Heliums and Frensley [fT2l using a matrix partitioning of the Hamiltonian that 
these boundary conditions are equivalent to an extra "memory" term - a time- 
dependent embedding or self-energy term - added to the Hamiltonian. This spatial 
partitioning, and the corresponding form of the embedding term, is appropriate for 
a Hamiltonian constructed with localized basis functions or spatially discretized 
wave-functions. 

An embedding method for solving the time-independent Schrodinger equa- 
tion in the region of interest was developed many years ago by the author [Tl3l[T4| . 
The operator which embeds region I - the embedding potential - is given by the 
surface inverse of a Green function for region II, evaluated over the boundary 
between the regions. Any convenient basis set can be used to expand the wave- 
functions (or Green function) in region I, and the method has been widely used 
in accurate surface calculations [15], and recently in photonics applications lfT6ll . 
The embedding potential is a generalized logarithmic derivative, giving the nor- 
mal derivative of the wave-function over the boundary of region I in terms of the 
amplitude - in mathematical terms it is a Dirichlet-to-Neumann map [fTTl . It is 
this form of embedding potential, rather than a tight-binding or discretized form, 
which we shall apply in this paper to the time-dependent problem. 

The time-dependent version of the Dirichlet-to-Neumann map has been stud- 
ied by Ehrhardt [Tl8l . and applied to embedding the time-dependent Schrodinger 
equation, discretized spatially as well as temporally - Ehrhardt paid particular at- 
tention to the stability and accuracy of the time-evolution. Moyer [fT9l has used 
this to study a range of one-dimensional problems including the scattering of elec- 
trons in model semiconductor structures. Recently Kurth et al. [20] have devel- 
oped a spatially discretized method to study quantum transport through a struc- 
ture between two leads, replacing the leads by time-dependent embedding self- 
energies; they also consider time- varying bias potentials in the leads. Boucke, 
Schmitz and Kull [[8]| have applied the time-dependent relationship between nor- 
mal derivative and amplitude to the one-dimensional problem of an oscillating 
potential of the form V(z,t) = — 1/ cosh 2 [^ + £o sin (cut)] acting on the wave- 
function 1 / cosh(^). (This wave-function is the bound state of the static potential; 
the oscillating potential is equivalent, by the Kramers-Henneberger transforma- 
tion [21], to the static potential in a spatially uniform time-varying electric field.) 
They solve this problem on a finite spatial grid, with the boundary relationship, 
for zero potential in the external region, applied to the end grid points. 
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Instead of using a finite-difference grid, Ermolaev et al. Il22l expand the wave- 
function in region I in a basis set. The matrix element of the kinetic energy op- 
erator evaluated over this finite region of space gives a normal derivative term at 
the boundary of region I, and this is evaluated using the Green function relation- 
ship with the boundary amplitude. They also consider the time-evolution of the 
bound state of the — 1/ cosh 2 (z) potential, but instead of transforming the spa- 
tially homogeneous, time-varying electric field using the Kramers-Henneberger 
transformation, they treat the problem directly, with a sinusoidally-varying field 
in region II. Like Ermolaev et al. we will also use a basis set for expanding the 
wave-function in region I, but we incorporate the amplitude-derivative relation- 
ship in a way analogous to the energy-dependent embedding method lfl~3l . 

We begin this paper with a derivation in section 2 of the time-dependent 
embedding method, based on our original method for embedding the stationary 
Schrodinger equation. In section 3 we derive the time-dependent embedding po- 
tential analytically for a constant external potential in region II, and demonstrate 
a numerical technique for embedding on to a constant electric field. We move on 
to model applications in sections 4 and 5. In section 4 we shall consider the same 
problem as Boucke et al. [%] - the oscillating potential — 1/ cosh 2 [z + £ sin(o;£)] 
acting on the 1 / cosh z bound state, using a very small region I enclosing this po- 
tential and embedding on to zero potential on either side. We shall see that the 
results are in excellent agreement with a finite-difference calculation extending a 
large distance on either side of the time-varying potential. In this section we shall 
next consider the excitation of a localized wave-function at a the surface of a free- 
electron metal by a time-dependent perturbation, with a static electric field on the 
vacuum side of the surface. We could use this in a model of photo-assisted field 
emission E31 . or even pump-probe experiments; here the calculation provides a 
test of the embedding potential for a uniform electric field. 

The problems we solve in section 4 are restricted to wave-functions confined 
initially to region I, though in the course of time-evolution they leak outside. This 
is a severe restriction for surface physics or other condensed matter applications, 
where the initial wave-function is usually an extended state. In section 5 we shall 
show how the time-evolution of such wave-functions can be calculated using em- 
bedding, with region I containing the time-dependent perturbing potential, but not 
necessarily the starting wave-function. The model calculation we shall describe 
in section 5 corresponds to the surface of a free-electron metal, with a sinusoidal 
time-dependent potential applied in the surface region (region I). With a small 
basis set, the time-evolution of the bulk wave-functions at the surface can be cal- 
culated very accurately. 

Atomic units are used throughout this paper, with e = h = m e = 1. The 
atomic unit of time = 2.418884 x 1(T 17 s, so that 1 fs = 41.34138 a.u. 
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2 Stationary and time-dependent embedding 

In this section we shall start from our previous results for embedding the station- 
ary Schrodinger equation lfT3l . and show how these can be transformed into time- 
dependent embedding. The embedding problem can be represented schematically 
by figure [TJ We wish to solve the Schrodinger equation - either time-independent 




Figure 1: Region I is embedded on to region II over surface S. 

or time-dependent - in the whole system, regions I + II: region I is treated ex- 
plicitly, and region II is replaced by an "embedding potential" at the interface S, 
added on to the Hamiltonian for region I. In a typical application to surfaces, re- 
gion I would be the surface layer or two of atoms together with the potential barrier 
region, and region II would be the vacuum on one side and the substrate crystal 
on the other [15J. The embedding potential ensures that the wave-functions eval- 
uated in region I match correctly in amplitude and derivative on to the appropriate 
solutions of the Schrodinger equation in region II. 

2.1 Stationary embedding 

The original embedding method is based on a variational principle for the energy, 
and we shall here outline the derivation lfT3l[T4l . The one-electron Hamiltonian, 
of which we wish to find the expectation value, is given by 

H = - l -V 2 + V{v) 1 (1) 
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where V is the one-electron potential. In region I we have an arbitrary trial 
function <p; this is extended through region II with the exact solution ip of the 
Schrodinger equation, evaluated at some trial energy e, which matches in ampli- 
tude on to 4> over the boundary S joining the regions (figured])- It is assumed that 
at an external boundary (the dashed line in figured]), the wave-functions satisfy a 
homogeneous boundary condition, typically going to zero. The expectation value 
of H is then given by 

E = Si d J^± + e Sn <hT± + \ Is dvsP (j| - g) 
fj dr<p*(f) + J u drip*ijj 

The first two terms in the numerator are the expectation values of the Hamiltonian 
in regions I and II. The third term, an integral over the boundary S, contains the 
difference in normal derivatives on either side of S (measured outwards from I) 
and comes from the kinetic energy operator acting on the kink in the trial function. 
We now use a result obtained by applying Green's theorem in region II, 

*(r,)=-i/ i *^(r„^; e )^ (3) 

- G is the Green function in II satisfying the zero normal derivative boundary 
condition on S. Taking the inverse of © gives the result which is central to the 
embedding method, 



dj>(r s ) 
dns 



-2 / dr f s G^(r s ,r' s ;e)iP(r f s ), (4) 
J s 



where Gq 1 is the surface inverse over S. Gq 1 is the embedding potential, and as 
we see from equation © it is a generalized logarithmic derivative. Gq 1 is the 
same as the mathematicians' Dirichlet-to-Neumann map IPT71 l24l . mapping the 
Dirichlet boundary condition (the amplitude specified over the boundary) on to 
the Neumann (the derivative specified). Using this result, and the fact that ip = <p 
over the boundary, we obtain 

J/^*-^- = -2 J s dr s f^sPiTgffi^TgJgieWg). (5) 

To complete the simplification of equation ©, we use a second result involving 
the embedding potential [TT3l . 

J n dr^ = ~J s dr s J s dv' s et>\v s ) dG ° 1{r d S ^ s]e U {v' s ). (6) 
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Substituting © and © into © gives us the embedding variational principle, 



E 



J, dr^H<p + \ J s dr s <P*£: + Js dr s J s dr^ G 



9Gn 



dnc 



fj dv(j)*(j) - f s dr s J s dr' 5 0*(r 5 



(7) 



- this gives E in terms of the trial function <fi defined in region I and on its bound- 
ary S. 

The wave-function <p which minimizes © satisfies the following equation in 
region I, including S, 

~V 2 + V(r)]o(r)4 



S(r- r s ) 



1 d(j) 



2 dn. 



+ (G \e) + (E - e 



9G 
de 



Ecj>{v). (8) 



The embedding potential evaluated at trial energy e plus the energy derivative term 
give Gq 1 at the required energy E, to first order in (E — e). The surface terms 
in square brackets vanish when has the correct normal derivative to match on 
to the solution of the Schrodinger equation in region II, so this is the correctly 
embedded solution of the Schrodinger equation. 

To find the solutions of © and ®, we expand <\> in terms of basis functions \i 
(here we assume that they are real and orthonormal), 



r . 



(9) 



Substituting into © and finding the stationary values of E gives the following 
embedded Schrodinger equation in matrix form, 



ijjdj 



Edi. 



(10) 



The first term comes from the Hamiltonian and the surface derivative term in the 
variational principle, 



Hi, 



drxi{r) 



while the second term comes from the embedding potential, 

dG^(v s ,r' s ;e) 



(11) 



o^s / dr' s Xi(r s ) 



G«\r s ,r' s ;e) + (E-e) 



de 



(12) 

When e lies in an energy continuum of region II, Gq 1 is complex, so we usually 
find the Green function of (H + E) rather than the eigenvectors. 
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2.2 Time-dependent embedding 

We shall now build the time-dependent formalism on the results of the last section. 
First we consider the the relationships between normal derivative and amplitude 
as functions of time, corresponding to © and © in the energy-dependent case. 
Taking the Fourier transform of © gives us 

$(r 8 ,t) = -\ f dv s f dt'G {v S: v' s -t-t')^^l (13) 
2 Js J~oo on s 

- we use a tilde to indicate functions of time. G satisfies the time-dependent 
Schrodinger equation in region II, 

(-^V 2 + V(t) - i£\ G (r, r'; t - t>) = S(r - r')S(t - f), (14) 

with the zero-derivative boundary condition on S, and because we take the re- 
tarded Green function we have, 

G (r,r';t-t') = 0, t < t' (15) 

- hence the upper limit of t in the integral in (TT3T) . This can be derived directly 
by applying Green's theorem to the inhomogeneous equation (TT4l) . and the corre- 
sponding homogeneous equation for ^(r, t) JH. 

Turning to the inverse relation giving the derivative in terms of the amplitude, 
we have to proceed differently, because the Fourier transform of Gq 1 does not 
converge - it is an increasing function of e. However, Ehrhardt [18J and Boucke 
et al. flU have shown that dip(t)/dns can be expressed in terms of the time- 
derivative of ij}(t) , and here we give a simple derivation. We re- write the integrand 
in © as 

J S 26 

- the Fourier transform of Gq 1 (e) / — ie converges, and the transform of —ie^p(e) 
is dil>{t)/dt. Defining Gq x (t) as the following Fourier transform, 

G \r s ,r' s ;t) = i- / + °° efeexp(-irf) °° ^ ^ ' } , (17) 
the Dirichlet-to-Neumann equation in time becomes 



drit 



-2 J dr s f dtGo^rs, r' 5 ; t - ' • (18) 
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We now turn to the embedding problem. The wave-function satisfying the 
time-dependent equation, at this stage with a time-independent potential, can be 
written in region I in terms of solutions of the embedded Schrodinger equation 



*) = £■ 



re 



-iEit 



(19) 



Equation ([8]) simplifies if the embedding potential is evaluated at the eigen-energy, 



V 2 + \/(r)U,(r) + 5(r-r 5 ) 



2dn< 



E; 



(20) 

and multiplying this equation by the coefficients in (fl9l) and summing over i gives 



-lv 2 + F(r))o(r./) 4 ,) { r-r s 



1 d(j) 



2dn s 

f dr' s £ OtGo 1 (r 5 , r' 5 ; ^^(^e"^* 
Using the same trick as in going from (fT6l) to (fTTT) . this becomes 



i- (21) 



V 2 + ^(r)U(r,t) + 5(r-r 5 ) 



1 90 



dr'o 



2dn s 



dt' 



(22) 



an equation which holds for r inside region I and on the boundary S. 

To show how the embedding terms in ((22l) work, we construct the solution of 
the time-dependent equation in region II, ip, with the inhomogeneous boundary 
condition that it matches in amplitude on to <fi over S, at all times up to t, 



4>(r s ,t) =0(r s ,f), t'<t. 
From (fTSl . the normal derivative on S is given by 



(23) 



dn,Q 



-2jjv s j dt'G \v s ,v' s -t-t')^^l. (24) 



But the right-hand side is the second term in the square brackets in (l22l) . with a 
factor of —2. For ((221) to be satisfied on S as well as inside region I, the two terms 
in the square brackets must cancel, forcing 0(r, t) to match in normal derivative 
on to the solution in II; as they already match in amplitude by construction, we 
have the correctly embedded solution of the time-dependent equation. Moreover, 
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the normal derivative term combined with —\V 2 give a Hermitian operator when 
integrating over region I. 

It is convenient to start off the time-evolution at t = 0, assuming that 

0(r 5 ,t) = O, t<0, (25) 



so we change the lower limit in the embedding term in (|22|) from — oo to t = 0. We 
also allow the potential in region I to be a function of time - we can presumably 
do this, as it does not affect the surface integral, originating from region II. This 
gives us the final form of the embedded time-dependent Schrodinger equation, 

V + V(r,t)Wt) + 5(r-r 5 ) 



2 



/ dr' s f' dt'G \r s ,r' s ;t-t 
Js Jo 



2dn s 
df 



>w (26) 



As in the time-independent case, this equation can be solved using a basis set 
expansion of 0(r, t), 

0(r,*) = $X*) Xi (r). (27) 



Substituting into (1261) we obtain the matrix form of the embedded Schrodinger 
equation, 

E (n«(tMt) + [dt% 3 {t - 1')^) = (28) 



j 

with the Hamiltonian matrix given by 



Hij{t) = i jf drVxi(r) • Vx,(r) + jdv X i{v)V{r,t) X ^), (29) 
and the embedding matrix by 

tij(t) = [ dr s [ dv' sX i{T S )G^{T S ,v' s -t)xjiy' s ). (30) 

J S J s 



The structure of (1261) and (1281) is the same as in the spatially discretized approach to 
the problem [fT2ll20l . with an embedding operator added on to the time-dependent 
Schrodinger equation. 

3 The time-dependent embedding kernel 

In this section we shall evaluate the embedding kernel Gq 1 (t) for constant and lin- 
ear potentials in region II. We work in one-dimension, but we shall see in section 
3.3 how the results can be applied to three-dimensional problems. 
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3.1 Embedding on to zero potential 

In the one-dimensional case we use (@J) to determine Gq 1 (e), and then substitute 
into the Fourier transform (fTTT) to find Gq 1 (t). With zero potential the wave- 
functions in region II satisfying outgoing boundary conditions for z > are 



^,e) = { eXP J~ 7 f' 7 = ^' £ n <0 (3D 
v ' \ exp(iKz), k = V2e, e > 0. 

Hence the embedding potential is 

G \e) = ^I/2 or -i^/2, (32) 

and substituting into (TTTT) the time-dependent embedding kernel is given by 

1 r+°° { ^— e < 

G \t) = — / deexp(-zet) (33) 

This integral can be done analytically [8J, and the result is 

We shall apply this embedding kernel in applications in sections 4 and 5. 



3.2 Embedding on to an electric field 

With other one-dimensional potentials it is necessary to carry out the Fourier trans- 
form in (TT7T) numerically, and we consider the linear potential corresponding to a 
constant electric field S. The wave-functions in region II satisfy the Schrodinger 
equation, 

ld 2 ip 

"2^-^ = ^' (35) 
which has Airy function solutions 11251 . 

M*) = Ai(-(2£)*(* + e/£)), 
and if) 2 (z) = Bi(-(2£)s(> + e/£)) . (36) 

Now we need the combination of ipi and Tp 2 corresponding to outgoing waves, and 
from the asymptotic behaviour of the Airy functions [25], this is given by 

■Mz) =Bi(-(2S)*(z + e/E))+iAi(-(2£) l s(z + e/E)) . (37) 
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So in this case the energy-dependent embedding potential at z = is 



Gr 



(2£)s 


Bi' 




+ iAi' 


f -(2£)h/Z) 


2 


Bi 




(2S)h/s) 


+ iAi 




(2£)h/£) 





(38) 



which we evaluate using the Airy function programs due to Gil et al. Il26l . For 
large positive or negative e, Gq l (e) has the y/e free-electron behaviour given by 
equation (|32|) . so once again we divide by e when finding the Fourier transform to 
obtain the time-dependent kernel. The free-electron behaviour at large e is what 
we would expect - the potential becomes irrelevant in this limit. 

In the numerical Fourier transform of Gq 1 (e)/e we must be careful about the 
singularity at e = 0, and we rewrite the integral as 



i r+oo 

Go {t) = 7T~ / deexp(-iet) 

2lX J-oo 



G H 



I 



2% J-oc e 2n J-oc e 

The first integral is well-behaved at e = 0, but in evaluating the second integral we 
have to be careful about the contour of integration around this point. As Gq 1 (t) is 
zero for t < 0, the singularities of the integrand must lie in the lower half-plane, 
so we replace ^ by where 77 is infinitesimal. Our Fourier transform then 
becomes 

^-1/^ i f + °° j / • .^o 1 ^) -GqHO) \ 0, t<0 

G " {t) = 2iL d ^-«t) U e ' + ^. i>a (40) 



To evaluate the integral in (1401) . our procedure is to include an exponential damp- 
ing term exp(— |e|/T), and discretize the integral with finite limits. This simple 
method of evaluating the Fourier transform works well. 

Results for the time-dependent embedding kernel with an electric field 8 = 2 
a.u. are shown in figure [2l using a coefficient Y = 500 a.u. in the damping term. 
We see that Gq is accurately zero for negative t. For t — > from above, the 
kernel tends to the free-electron value (l34l) in accordance with our intuition: at 
short times (or high energies) an electron cannot tell whether it is in an electric 
field or not. For larger t, Gq 1 (t) rapidly approaches the zero-frequency embed- 
ding potential, Gq 1 (0), which is the contribution from the pole in the frequency 
integral (|39| ); the larger the field, the more rapidly it approaches a constant value. 
This form of Gq 1 (t) makes it very easy to evaluate the long-time contribution to 
the embedded Schrodinger equation. 

We test the accuracy of the embedding kernel by calculating the time-evolution 
of a wave-function ip(z, t) over the extended system of regions I and II, and then 
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Figure 2: G 1 (t) for the Schrodinger equation with electric field £ 
line, 1ZeG 1 ; dashed line, TrnGg 1 . 



2 a.u. Solid 



use (|18l) to compare the Dirichlet-to-Neumann result with the directly-calculated 
derivative at the boundary of region II. In this test the potential itself is time- 
independent, with an infinite barrier at z = 0, zero potential between z = 
and z', and in region II beyond z' an electric field corresponding to the potential 
—£(z — z'). The wave-function at t = is taken to be a sum of Gaussians 
vanishing at z = 0, 



ip(z, 0) = exp 



Z() 



w 



exp 



z + z 



w 



z>0, 



(41) 



and the value of z and the width w are chosen so that the initial value of the 
wave-function at z', where the electric field starts, is negligible. To solve the time- 
dependent Schrodinger equation for ip(z, t) we discretize in both space and time, 
evolving forward in time steps St using the norm-conserving Crank-Nicolson 
method BH, 



ip(z,t + 6t) 



1 + 



iStH 



' iStH\ ~, v 
l- — U(z,t), 



(42) 



where H is the finite-difference Hamiltonian matrix QUI . 

For this test we take an electric field 8 = 2 a.u. beyond z' = 3 a.u., with the 
constants in the Gaussian wave-function z = 0.5 a.u., w = 1 a.u. The intervals 
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Figure 3: Electric field test: dip /dz at z' evaluated from the Dirichlet-to-Neumann 
relation compared with direct calculation from the numerical wave-function. D- 
to-N results are the solid line (real part) and the dashed line (imaginary part); 
direct results are the dots, hardly distinguishable except close to t — 0. 

in the spatial and time discretization are 5z = 0.0005 a.u. and 5t = 0.005 a.u., 
and the spatial range extends to z = 40 a.u. - enough to eliminate edge effects 
over the time-range we consider. We then evaluate dip/dz at z' directly from 
the finite difference results, using first order differences, and compare this with 
the derivative evaluated using (fT8l) with the embedding kernel shown in figure [2l 
and first-order differences for dip/dt. The comparison is shown in figure|3l the 
values of dip/dz from the Dirichlet-to-Neumann relationship (solid and dashed 
lines) can hardly be distinguished from the results for the derivative calculated 
directly (dots). Only at very short times, and because of some noise in the direct 
results, can they be distinguished. This shows that our numerical evaluation of the 
embedding kernel for the electric field is accurate - the method can presumably 
be extended to other potentials. 
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3.3 Shifting the potential 

If we know the embedding kernel for a one-dimensional potential V(r) in region 
II, we can immediately find an expression analogous to the Dirichlet-to-Neumann 
relationship (fTST) for the potential shifted by a constant V . The time-dependent 
Schrodinger equation satisfied by ip(r, t) in region II is given by 



X 



--V 2 + V(v)+V Q ^(r,t)=i-^, (43) 

and making the substitution 

#(r,t)=exp(-iV t)%(r,t), (44) 
^ satisfies the unshifted Schrodinger equation 

±V 2 + V(r)^(r,t)=i^. (45) 

As /dns is related to d^f/dt by the original embedding kernel for V(r), the 
corresponding relationship for ip with the shifted potential is given by 

— = -2exp(-tV t) 

dn s 

j s dr s J* dt'G \r s , r' s ; t - t')^[exp(iV t')i>{v' s , 0]- (46) 

This result will be particularly useful when we come to deal with the potential 
step at a surface. Moreover, a one-dimensional potential has often been used in 
surface calculations [[27], and then the solutions of the Schrodinger equation have 
the form 

4>(r, t) = exp(2K ■ R)j> K (z, t), (47) 

where K is the Bloch wave-vector parallel to the surface, and R is the surface- 
parallel component of r. Then ipK satisfies the one-dimensional Schrodinger 
equation shifted by K 2 /2, and the Dirichlet-to-Neumann relationship becomes 

d^Kiz.t) n ( iKH\ ft u A 9 [ (iKH'\ 7 , 

= -2 ex P j I dt>G Q \t - t>)- [ex P ^_ j ^ t>) 

(48) 

where Gq 1 is the embedding kernel for the one-dimensional potential. The right- 
hand side of (l48l provides the embedding potential for states with Bloch wavevec- 
tor K in this form of potential. 

Another approximation for surfaces is to assume the full three-dimensional 
potential for the surface region itself (region I), and a one-dimensional potential 
for the bulk substrate and the vacuum (which together constitute region II) 11281 . 
Again (l48~l) can be used to construct the embedding potential for each Fourier 
component of the wave-function in region I. We shall explore this in a later paper. 
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4 Embedding 

4.1 Model atomic problem 









-d 







+d +E 






Figure 4: Oscillating model atomic potential. Region I, treated explicitly, lies 
within z = ±d. The basis functions for expanding the wave-function in region I, 
are defined in terms of z = ±D. 

To test our time-dependent embedding method with the kernels derived in 
section 3, in this section we calculate the time-evolution of states which are ini- 
tially localized in region I. We start off with the time-evolution of the normalized 
bound state wave-function ^p(z) = l/y/2 cosh(,z) of the one-dimensional poten- 
tial V{z) = — 1/ cosh 2 (z) in a time- varying electric field S = So sin u>t. Ermolaev 



et al. [|22| solve this problem directly, but we follow Boucke et al. [8], who use 
the Kramers-Henneberger transformation 112711 to convert this into the problem of 
the wave-function evolving in the oscillating potential, 



V(z,t) = -l/cosh 2 [;z + £ sin(u;t)] 



(49) 



with zero potential beyond (figureHl). Here £ is the classical amplitude of oscilla- 
tion in the electric field, £ = S /u 2 . We solve the embedded Schrodinger equa- 
tion (I26ll28l) in region I, defined as \z\ < d, replacing the regions with \z\ > d by 
the zero-potential embedding kernel. The time-dependent wave-function in region 
I, 4>(z, t), is expanded in a basis set given by 



cos 2 !^, m even 



sm 



inn z 
2D ' 



m odd 



(50) 
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where D lies beyond d (figure HI) to give flexibility in amplitude and derivative 
at the boundary of region I. We orthogonalize and normalize the basis functions 
within region I by diagonalizing the overlap matrix S mn , 

S mn = dzE m (z)E n (z). (51) 



Defining a? as the j'th eigenvector of S, with eigenvalue Sj, the j'th orthonormal- 
ized basis function is given by 

Xj(z) = -^Y, a] m^m{z). (52) 

V s j m 

If overcompleteness is a problem, this will show itself as a very small value of Sj, 
and the corresponding basis function can be dropped. 

The time-dependent matrix equation (|28l becomes, in this one-dimensional 
case with embedding at both ends of the range, 



,.d<j>(-d, f) 



+ Xi(-d) f Q dt'G,\t - t, m , 

+Xi(d) J o dtG (t-t) — — — = l -fai ( 53 ) 
where the Hamiltonian matrix element is given by 

H tl (t) = j^dz {l^^ +Xi { z )V{z,t)xj{z)^ , (54) 

and the embedding kernel Gq 1 by (|34|) . 

We turn to the numerical time-integration of (1531) . which we write as 

^ + iHa = -iT, (55) 
dt 

where T is the vector representing the embedding terms in (1531) . 



Ti = Xi{-d) j dtG l {t-t ) — + Xi{d) j dt G r (t - t ) — ^— . 

(56) 

We can improve on the first-order integration scheme, 

a(t + 5t) = [1 + idtHit)}- 1 ^) - idtT(t)], (57) 
by expanding the time-evolution operator to second order in St, giving 

CSt 2 H(t) 2 \ 1 
1 + iStH{t) I [a(t) - iStT(t)]. (58) 
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Although it is not consistently second order in St, this stable scheme proves more 
accurate than (1571) in our tests. In evaluating the integrals in ((56l) we use 



d<j>{t) (59) 



at st 

rather than the more accurate formula 

d<t>(t) <j>(t + St) - 4>(t - St) 



(60) 



dt 2St 

because (l60l) cannot be applied at the upper limit of the integral - we do not 
yet know <p(t + St). To do the integration itself we subtract off the (t — t')^ 1 / 2 
singularity at t' = t, and then use the trapezium rule for the remaining well- 
behaved integral. 

We calculate (p(z, t) in the oscillating potential (|49l ), with an amplitude of 
oscillation £ = 2.5 a.u. and frequency u> = 0.2 a.u., corresponding to an electric 
field £ = 0.1 a.u. (These are the values used by Boucke et al. [8] - as the 
bound state energy of the starting wave-function is —0.5 a.u., ionization is due to 
multiphoton processes.) Region I is taken with d = 10 a.u., and the basis functions 
are defined with D = 13 a.u.; the results presented below are for 25 and 40 basis 
functions. An interval St = 0.01 a.u. is used in the time integration (|58l) . As 
a benchmark we compare the embedding results with the wave-function ip(z,t) 
calculated over an extended range, using finite differences and Crank-Nicolson 
time-integration (l42l) (throughout this section we use ip(z, t) to indicate the wave- 
function through space extended beyond region I). The extended wave-function is 
calculated over the range —400 < z < 400 a.u., with a spatial interval Sz = 0.004 
a.u., and a time-interval St = 0.02 a.u. 

Figures [51 [6] and |7J show the comparison between the embedded and finite 
difference wave-functions in region I, for times t = 80 a.u., 320 a.u., and 400 a.u. 
(the period of the oscillating potential is 31.42 a.u.) and we see that the agreement 
is very good. At the scale of figures [5ta)-|6ta), the magnitude of the embedded 
wave-function \<p\, calculated with 40 basis functions, is indistinguishable from 
the magnitude of the extended wave-function From figures Ob)-[6tb), we see 
that the error in \<p\ is about 5 x 10~ 4 a.u., with a slight decrease in accuracy 
with increasing time. The results with 25 basis functions are less accurate - they 
are just distinguishable from the extended results on the scale of figures [5fa) and 
[6ta). At t = 400 a.u. shows slight oscillations about the embedded \<p\ with 
40 basis functions (figure [7J a)). The reason for this is that the extended wave- 
function has been reflected from the limits of its range at z = ±400 a.u., with 
resulting interference, as we can see from figure iTjb). Our embedding results are 
less accurate than those of Boucke et al. [H, who achieve a relative accuracy 
in their embedded wave-function of about 5 x 10~ 5 ; however, they apply their 
embedding as a boundary condition on a much larger finite-difference calculation. 
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Figure 5: Oscillating potential: magnitude of wave-function at t = 80 a.u., 2.55 
periods. |0| calculated using embedding: solid line, 40 basis functions; dashed 
line, 25 basis functions. calculated over extended space using finite differ- 
ences: short-dashed line. 

(a) Plotted over embedding region between z — ±10 a.u. |0| with 40 basis func- 
tions and |^| are indistinguishable on this scale. 

(b) Plotted around the embedding point at z = 10 a.u. 
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Figure 6: Oscillating potential: magnitude of wave-function at t = 320 a.u., 10.19 
periods. |0| calculated using embedding: solid line, 40 basis functions; dashed 
line, 25 basis functions. calculated over extended space using finite differ- 
ences: short-dashed line. 

(a) Plotted over embedding region between z — ±10 a.u. |0| with 40 basis func- 
tions and |^| are indistinguishable on this scale. 

(b) Plotted around the embedding point at z = 10 a.u. 
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Figure 7: Oscillating potential: magnitude of wave-function at t = 400 a.u., 12.73 
periods. 

(a) \<f>\ calculated using embedding: solid line, 40 basis functions; dashed line, 25 
basis functions. calculated overextended space using finite differences: short- 
dashed line. Plotted over embedding region between z = ±10 a.u. 

(b) 1^1 plotted over extended region between z = ±400 a.u., showing reflection 
of the wave-function at the boundaries of the finite-difference calculation, and 
subsequent interference. 
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4.2 Surface excitation in a field 



VJI 





Figure 8: Model atomic potential near a jellium surface, with applied electric field 
in the vacuum. Region I, treated explicitly, lies within embedding boundaries at 

zi and z r . 

In this section we shall follow the time-evolution of an electron wave-function, 
initially in the l/y / 2cosh(2; — z') bound state of a —1/ cosh 2 (;z — z') potential 
near a surface, at which there is a constant applied electric field £ (figure [8]). The 
surface potential step V st is broadened, so the static potential felt by the electron is 
given by 



where 9(z) is the step function. The surface parameters we use are V st = 0.5 a.u., 
£ = 0.5 a.u., and we position the "atomic" potential at z' = —4 a.u. The electric 
field is £ = 0.2 a.u., approximately 10 11 Vm _1 - this is much larger than the fields 
used in field emission but it could represent the field of an IR laser in a quasi-static 
approximation. To excite the bound state electron, an additional time-dependent 
potential is applied, of the form 



with a = 1 a.u., ( = 2 a.u., and frequency u = 1 a.u. This is turned on at 
t = 0, and we follow the subsequent time-evolution of the bound-state electron 
wave-function. 




(61) 




(62) 
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Region I lies in the surface region between the embedding boundaries at zi 
and z r (figure [8]). We take Z\ = — 14 a.u., and z r = 6 a.u., so the embedding 
region extends ±10 a.u. on either side of the atomic potential. The Hamiltonian 
for region I is embedded on to the zero-potential embedding kernel (1341) at zi ; at z r 
it is embedded on to the Airy function kernel evaluated numerically (section 3.2), 
shifted by the constant potential V (figure [8]), using the shift formula (|46l) . We 
use the basis functions given by (l50l) to expand the time-dependent wave-function 
in region I - these are centred on the atomic potential with D = 13 a.u. 

Results for <fi(z, t) at t — 50 a.u. are shown in figure EJa), calculated with 
25 basis functions, compared with a finite-difference calculation for ip(z,t) cal- 
culated over an extended range, —400 < z < 400 a.u. (figure Hfb)). We see 
that agreement is very satisfactory, though not quite as good as in section 4.1 
where we used an analytic embedding kernel. In this calculation we need to take 
St = 0.0025 a.u. in the time integration, compared with 0.01 a.u. in section 4.1. 
This is probably because of the difference in the time-dependent perturbation. 

We also calculate the current crossing the embedding boundaries at zi and z r , 
using the expression 

j(z,t) =Zm($*^Y (63) 

with the normal derivative determined from the embedding formula (fT8l) . Taking 
the derivative outward from region I, a positive current indicates charge leaving 
the region. Our results are shown in figure [ToTa) for the current crossing each 
boundary as a function of time - here we have something physical, which could 
in principle be compared with experiment. There is excellent agreement with the 
current calculated from the finite difference results, as we see from figure [TOlb). 
giving the current in the time range where there is the biggest difference between 
the two methods. 



5 Time-evolution of extended states 

The formalism developed up to now assumes that the wave-function 0(r, t), whose 
time evolution we study in region I, has zero amplitude for t < at the embed- 
ding surface and in region II. But in condensed matter applications we are usually 
interested in exciting bulk states to which this condition does not apply, and to 
study their time-evolution we have to extend the formalism. 

We start with a wave-function \I/(r) which is an eigenstate with energy E of 
the time-independent Hamiltonian H , extending through regions I and II. For 
times t > a time-dependent perturbing potential 5V(r, t) is applied - confined 
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Figure 9: Excitation in a field, £ = 0.2 a.u.: magnitude of wave-function at t = 50 
a.u. 

(a) \<j)\ calculated using embedding with 25 basis functions, solid line; cal- 
culated over extended space using finite differences, dashed line. Plotted over 
embedding region between z = — 14 a.u. and 6 a.u. 

(b) 1*01 plotted over extended region. 
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Figure 10: Excitation in a field, S = 0.2 a.u.: current across boundaries of region 
I as a function of time. 

(a) j calculated with embedding: solid line, across left-hand boundary at z\\ 
dashed line, across right-hand boundary at z r . 

(b) Comparison of j at z r using embedding (solid line), with finite differences 
(dashed line). 
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to region I - and the wave-function is subsequently given throughout space by 

$(r, t) = *(r) exp(-iEt) + fj{r, t). (64) 
Substituting into the time-dependent Schrodinger equation gives 

[H + 5V(t)]fj(t) + 5V(t)Vexp(-iEt) = (65) 

- fj(r, t) satisfies the time-dependent Schrodinger equation with an additional in- 
homogeneous term. In region II, where 5V = 0, this term vanishes and fj satisfies 
the original Schrodinger equation, 

H fj(t) = i^t, (66) 

so fj and its normal derivative across S are related by the Dirichlet-to-Neumann 
result (fT8l) . 

dfj(r s ,t) r t l , + ,n-x, ^dfj(r' s ,t' 



dn<. 



J s dr s J o dt'G^vs^t-t') V Y t ; (67) 



This means that we can write a time-dependent Schrodinger equation for fj (r, t) 
analogous to ((261) . with the extra inhomogeneous term, 



if. (68) 
at 



+8(r - r s 



X - V 2 + Vb(r) + 5V(r, tj) fj(r, t) + SV(r, (r) exp(-iEt) 



Solving this equation within region I gives us the change in bulk wave-function. 
Note that in this section we use ^(r, t) for the change in wave-function with r in 
region I, evaluated by embedding (analogous to <p m previous sections), as well as 
the change in wave-function extended through all space (analogous to iff). 

As a model problem we consider the jellium surface with a smeared-out po- 
tential step, 

V4W = V. ( 1 + t 7"^ ) , (69) 

using the same parameters as in section 4.2, V st = 0.5 a.u., £ = 0.5 a.u. Region 
I lies within embedding boundaries at zi = —10 a.u. and z r = 10 a.u., beyond 
which the potential is taken as constant. The bulk continuum wave-function is 
found numerically using Numerov's method [30], matching W(z) on to the asymp- 
totic solutions of the time-independent Schrodinger equation 

^ (sin(^ + 0) z< Zl 
I aexp(— 72;), z > z r 
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with 



k = V2E, 1 = ^{2{V A -E). 



(71) 



We use the same time-dependent perturbing potential as in section 4.2 (l62l) . but 
in this case centred at the surface step. Region I is embedded at Z\ on to the free- 
electron embedding potential given by (1341) . and to the right on to this embedding 
potential shifted by the potential step V st , following (l46l) . 

The expansion coefficients in the basis set expansion of fj(z, t) in region I 
satisfy (l53l) and (l55l) . with an extra term b{ added on to r, ; (l56l) given by 



As in the surface electric field calculation in section 4.2, we take 5t = 0.0025 a.u. 
in the time-integration for accurate results. 

We take a bulk state with energy E = 0.3 a.u. and apply the time-dependent 
surface perturbation with frequency uj = 0.5 a.u., starting at t = 0. The results 
for the modulus of the change in wave-function at t = 120 a.u., calculated with 
25 basis functions, are shown in figure ITTT a), compared with results from a finite- 
difference calculation taken over the extended range —1000 < z < 1000 a.u. 
Figure \TT\b) shows the finite difference results over part of the extended range, 
and we see edge effects, propagating from the left-hand boundary and reaching 
z = —100 a.u. With the left-hand boundary in the finite difference calculation at 
— 1000 a.u., the results in the surface region become unusable beyond t = 120 a.u., 
a problem which of course does not affect embedding with its correct treatment 
of boundary conditions. We see from figure 1TTT a) that the embedding results are 
accurate, even with this relatively small basis set. 

It is interesting to calculate the current in this case, not only as a sensitive test 
of the accuracy of the calculation, but also to illustrate the physics. Of course 
we must use the full wave-function given by (l64l) in the expression (l63l) for the 
current. Checking the accuracy first, we see from figure [TUT)) that embedding 
works well. It is figure |T2t a) which has physical content, and we see that after a 
short period of transient behaviour, the current entering the surface from the bulk 
(a negative current at zj) and leaving the surface into vacuum (a positive current 
at z r ) both settle down to steady behaviour. In fact the currents balance out on 
average, as is shown by figure [131 giving the norm of ip(z,t) in the embedding 
region. We see that the charge in the surface region oscillates about a constant 
value, after a sudden loss of charge at about t = 10 a.u. 

The results shown in figure [T2T a) are in some ways surprising - we are solving 
the Schrodinger equation in the surface region with an embedding potential based 
on an outgoing Green function, and yet we are able to describe the current entering 
the surface from the bulk. This goes to show the power of Green function-based 
methods. 




(72) 
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Figure 11: Excitation of a bulk state: magnitude of change in wave-function at 
t = 120 a.u. 

(a) 1 77 1 calculated using embedding with 25 basis functions, solid line; \fj\ cal- 
culated over extended space using finite differences, dashed line. Plotted over 
embedding region between z = —10 a.u. and 10 a.u. 

(b) I fj I from finite differences plotted over extended region. 
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Figure 12: Excitation of a bulk state: current across boundaries of region I as a 
function of time. Positive current corresponds to charge leaving region I. 

(a) j calculated with embedding: solid line, across left-hand boundary at z\\ 
dashed line, across right-hand boundary at z r . 

(b) Comparison of j at z r using embedding (solid line), with finite differences 
(dashed line). 
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Figure 13: Norm of wave-function in the embedding region as a function of time. 

6 Conclusions 

This embedding method provides the correct boundary conditions for solving the 
time-dependent Schrodinger equation in a limited region of space, region I, auto- 
matically matching the solution on to the time-evolving wave-function in the rest 
of the system, region II. Once we have found the embedding kernel for region II, 
all that we have to do is to add this on to the Hamiltonian for region I and time- 
integrate the Schrodinger equation in this region. As we solve the time-dependent 
Schrodinger equation using the relatively small basis set needed to describe region 
I, this embedding method is very economical. In the examples in this paper we 
use a plane wave basis to expand the the wave-function in region I, but any other 
basis set (see, for example ref. PTlO could be used. 

The next stage in this project is to improve the numerical time-integration 
scheme, and then apply it to more realistic surface models to study surface electron 
excitation. 
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